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Abstract. Multiscale dynamics are frequently present in real-world processes, 
such as the atmosphere-ocean and climate science. Because of time scale sepa- 
ration between a small set of slowly evolving variables and much larger set of 
rapidly changing variables, direct numerical simulations of such systems are dif- 
ficult to carry out due to many dynamical variables and the need for an extremely 
small time discretization step to resolve fast dynamics. One of the common reme- 
dies for that is to approximate a multiscale dynamical systems by a closed approx- 
imate model for slow variables alone, which reduces the total effective dimension 
of the phase space of dynamics, as well as allows for a longer time discretiza- 
tion step. Recently we developed a new method for constructing a deterministic 
reduced model of multiscale dynamics where coupling terms were parameter- 
ized via the Fluctuation-Dissipation theorem. In this work we further improve 
this previously developed method for deterministic reduced models of multiscale 
dynamics by introducing a new method for parameterizing slow-fast interactions 
through additive stochastic noise in a systematic fashion. For the two-scale Lorenz 
96 system with linear coupling, we demonstrate that the new method is able to 
recover additional features of multiscale dynamics in a stochastically forced re- 
duced model, which the previously developed deterministic method could not 
reproduce. 



1. Introduction 

Multiscale dynamics are common in applications of contemporary science, such 
as geophysical science and climate change prediction |12[|13[|T9ll21[|24ll28[|35ll . Mul- 
tiscale dynamics are typically characterized by the time and space scale separation 
of patterns of motion, with fewer slowly evolving variables and much larger set of 
faster evolving variables. This time-space scale separation often made direct nu- 
merical computation of the dynamics on the slow time scale quite difficult in real- 
world applications, which led to the development of multiscale computational 
methods 111511161 . These methods make use of the averaging formalism ||36ll4Tll42l 
to allow for large time discretization steps for the computation of the slow part 
of the dynamics. However, in very large systems with many fast variables even 
these methods are computationally expensive. 
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As a different alternative to direct numerical simulation of the complete multi- 
scale model with all variables, it has long been recognized that, if a closed simpli- 
fied model for the slow variables alone was available, one could use this closed 
slow-variable model instead to simulate the statistics of the slow variables. In 
order to derive such a reduced model, one usually represents the process as a 
general two-scale dynamical system of the form 

(1.D £ = f t =\^,y), 

where x G R N - Y and y G IR N y are the state vectors, and F and G are nonlinear 
differentiable functions. The scaling parameter £ <C 1 is used to separate the 
time scales in into slow (that is, x), and fast (that is, y). The scale separa- 



tion parameter e above is introduced for convenience of presentation, because, as 
shown below, the method developed here does not require its explicit presence to 
function. Besides, often real-world processes do not have any distinct time-scale 
separation parameters, and the fast and slow variables are known empirically 
from observations. 

Under the assumption of "infinitely fast" t/-variables, one writes the reduced 
averaged system for slow variables alone as 

dx - - I' 

(1.2) — = F(x), F(x) = F(x,y)dMy), 

where is the invariant distribution measure of the uncoupled fast dynamics 
where x is a fixed parameter: 

(1.3) ^ = G(X,z). 

Under the assumption of ergodic jif, the measure integral in (|1.2|) can be replaced 
with the time average along a single trajectory Zx{t) of (|1.3|) for specified parame- 
ter x: 

(1.4) F(x) = lim I [ T F(x,Zi{t))dt. 

T— >oo 1 Jo 

From (|1.2|) - (|1.4|) , it is clear that the computation of the averaged function F(x) is 
not a simple task; in fact, it is rarely available explicitly (except for some special 
cases where either the invariant distribution measure or the solution z*(f) are 
known as explicit formulas). A good example of the system where the invariant 
distribution measure is known explicitly is the Ornstein-Uhlenbeck process l|40ll ; 
the invariant distribution measure of the Ornstein-Uhlenbeck process is a Gauss- 
ian distribution with explicitly known mean state and covariance matrix. Because 
of this convenient property, the Ornstein-Uhlenbeck process is popular in the area 
of stochastic modeling of geophysical and other real-world processes. 
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The "upgrade" from the deterministic reduced model in (|1.2|) is a stochastic 
model of the form 

(1.5) dx = F(x)dt + o-(x)dW tf 

where Wf is a Wiener process, and c(x) is a matrix. The stochastic reduced 
model of the form (|1.5[) is generally considered to be more advantageous to the 
deterministic reduced model in (|1.2[) , because the random noise can be used to pa- 
rameterize "noise-like" influence from fast variables onto slow dynamics, which 
is entirely lacking in the deterministic reduced model in (|1.2|) . Sometimes, de- 
terministic reduced models of atmospheric processes fail to adequately represent 
atmospheric variability; it is thought that a significant portion of atmospheric 
variability is essentially noise-like, and fitting it with deterministic forcing terms 
does not seem to be a good way to produce an adequate approximation. Ad- 
ditionally, the deterministic slow dynamics in (|1.2|) tend to produce the attractor 
of the slow dynamics with significantly lower dimension than that for the full 
two-scale dynamics in especially for weakly chaotic slow variables. At the 
same time, if the stochastic forcing in (|1.5[) is strong enough, it should "diffuse" 
the invariant manifolds of (|1.2[) , thus inflating the dimension of the set containing 
limiting dynamics. 

Nonetheless, in many applications, the averaged dynamics of the deterministic 
type (|1.2|) or stochastic type (|1.5[) are not available explicitly. As a result, nu- 
merous approximate closure schemes were developed for multiscale dynamical 
systems | [Tll[T6ll30l - l33l , which are all based on the averaging principle over the 
fast variables j|36ll41ll42~|] . Some of the methods (such as those in H301-33I) replace 
the fast nonlinear dynamics with suitable stochastic processes [44], discontinuous 
Markov jump processes 1122 II , or conditional Markov chains fl4f, while others Ifl6ll 
provide direct closure by suitable tabulation and curve fitting. Reduced stochas- 
tic dynamics were used to model global circulation patterns Illl[|18il34[|43[|45l , and 
large-scale features of tropical convection II23U291 . However, it seems that all these 
approaches require either extensive computations to produce a closed model (for 
example, [14111611 require multiple simulations of fast variables alone with different 
fixed states of slow variables), or somewhat ad hoc determination of closure coef- 
ficients by matching areas under the time correlation functions ||30] - |33| . Another 
interesting method was recently developed in MIDI , however, again, somewhat 
ad hoc approach was used to compute the closure (namely, slow variables of a 
multiscale system were treated as if their dynamics were generated by a running 
average of an Ornstein-Uhlenbeck process). 

In the two recent works |3H6l the author developed a relatively simple and 
straightforward method of constructing the deterministic reduced model for slow 
variables of a multiscale model with nonlinear and multiplicative coupling, which 
required only a single computation of certain statistics of the fast dynamics with 
a fixed state of the slow variables, located in the region where the slow dynamics 
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usually evolve. The method was based on the first-order Taylor expansion of the 
averaged coupling term with respect to the slow variables, which was computed 
using the Fluctuation-Dissipation theorem ||T]-l3ll7l-l9l l27ll38| . It was demonstrated 
through the computations with the appropriately rescaled two-scale Lorenz 96 
model ||25ll26l that, with nonlinear and multiplicative coupling in both slow and 
fast variables, the developed reduced model produced good approximation to 
the statistics of the full two-scale Lorenz 96 model. Among the advantages of 
the developed method were its simplicity and explicit formulation. Additionally, 
existing zero-order models of this kind for the Earth's atmosphere (such as the 
T21 barotropic model ||9HT7l|39l) can be retrofitted with the new deterministic 
correction term emerging from the theory in |4J. 

In the current work, we introduce a method for computing a consistent approx- 
imation of the form (|1.5|) to the multiscale dynamics in The method is aimed 
at stochastic parameterization of general complex nonlinear multiscale dynamics 
with many variables, and has essentially the same implementation restrictions as 
the method for deterministic reduced models we developed previously in BH6]]. 
We test the new method on the two-scale Lorenz 96 model, where only the linear 
part of the coupling is enabled for the simplicity of presentation. We demonstrate 
through direct numerical simulations that the new method generally improves 
properties of statistics of the reduced model, and, in particular, the injection of 
random noise into a deterministic reduced model can make it both more or less 
chaotic and mixing, depending on the difference between the dynamical regimes 
of the reduced and full two-scale models. 

The manuscript is organized as follows. In Section |2] we present the general 
description of the new method, following the homogenization theory of II37L In 
Section |3] we lay out the step-by-step computational implementation of the new 
method for a general two-scale dynamical system with linear coupling between 
the slow and fast variables, which does not contain any explicit time-scale sepa- 
ration parameters. In Section H] we test the new method on the two-scale Lorenz 
96 model in a range of dynamical regimes with varying chaos, mixing, and time 
scale separation between the slow and fast variables. Section [5] summarizes the 
results of this work. 



2. General description of the method 



Here we introduce a new method for the stochastic correction of the form (|1.5|) 
for the deterministic reduced model (|1.2|) . To derive the method, we use the theo- 
retical framework similar to that applied to homogenization problems in 11371 . To 
simplify presentation, we assume that the deterministic averaged function F(x) 
from (|1.2|) is already available, either as an explicit formula, or as an approxi- 
mation we developed previously in BH6J. Below, x is used to denote the state 
of the slow variables from the multiscale system while x. denotes the state 
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of the slow variables from the deterministic reduced model (|1.2|) . The difference 
between x and x is denoted as q = x — x. Then, one can rewrite the multiscale 



system in in the new variables as 

dx 
~dt 



(2.1a) 



-^- = F(x) 



(2.1b) 



(2.1c) 



^L = F(x + q,y)-F( 

dy 1 _ . . 
— = -G{* + q,y). 



What we see above is that the first equation is already a closed system from (|1.2[) , 
and given a suitable approximation for F from flUm, it can be solved on its own. 
Thus, x(t) can be treated as a given function of time. This, in effect, leaves q and 



y as the unknown variables, and the first equation in ((2JJ can be dropped. Now, 
the idea is to apply the averaging formalism to q, obtaining q as the next order 
correction to x. However, the straightforward application of what was done in 



(|1.2|) to (|2.1|) leads to q being identically zero for all times as long as its starting 
value is zero. In this situation, the reduced model for q should be derived as 
the Ito diffusion process of the corresponding backward Kolmogorov equation 
restricted to slow time scale (see [37] for details), under the condition that q (and, 



therefore, dql dt) in (|2.1[) is 0(e). To do that, first we bring the time scale of q in 
(|2.1|) to 0(1) by rescaling the time t as x = et. For the rescaled time x, from (|2.1[) 
we obtain 



(2.2a) 



(2.2b) 



^L = l[F(x + q,y)-F(x)], 

dy 1 N 
— = -G(x + q,y). 



Above, the evolution equation for x is no longer needed, as x(x) is a given 
function of rescaled time T. Now we write the backward Kolmogorov equa- 
tion for (|2.2|) . For that, let v(x, x',x, q, y), x' > x, be the value of a test function 
h(q(x'),y(x')), given q(x) = q, and y(x) = y. Then, v(x,x',x,q,y) obeys the 
following backward Kolmogorov equation (for reference, see, for example, [20J): 
(2.3) 



dv{x, x',x,q, y) 
dx 



1 1 - 

-^G(x + q r y) ■ V y + - [F(x + q,y) - F(x)] ■ V< 



v(x, x',x,q,y). 



Below, we drop the r'-dependence from v as it is of no consequence to what is 
presented. Observe that the terms in the backward Kolmogorov equation above 
are multiplied by different powers of e, which leads to the perturbation expan- 
sion of the solution v(x,x,q,y) in powers of e, and derivation of the backward 
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Kolmogorov equation for the term which has the lowest power of e in the closed 
form. Then, its corresponding Ito diffusion process will be the evolution equation 
for dql dt, as long as q is small enough. 

The expansion of v(t, x, q, y) in powers of e is 

(2.4) v = v + evi + e 2 v 2 + .... 

Plugging the expansion above back into the Kolmogorov equation in (|2.3)> and 
collecting the terms with matching powers of e, we obtain the following relations 
for each power of e: 

(2.5a) G(x + q,y) ■ V y v = for e~ 2 , 

(2.5b) G(x + q,y) -VyVi = -[F(x + q,y) - F(x)] -VqVo for e -1 , 



(2.5c) ^ = G(x + q,y)-V y v 2 +[F(x + q,y)-F(x)]-V q v 1 for e° 



Below we consider each of the relations above separately. 

• Order e~ 2 . From the relation in (|2.5a[) we determine that Vq(t,x, q,y) = 
Vo(T,x,q), that is, vq does not depend on y. 

• Order e _1 . Here we use the relation in (|2.5b|) to express V\ in terms of Vq. 
We denote the flow, generated by (|1.3|) , by (p s x , so that cp^y is the solution 
of (|1.3|) forward in time s with the initial condition y, with the obvious 
identity 

(2-6) ds^+^ = G (x + c l> fx+qV) ■ 

Now, consider the integral 



(2.7) 



u(s,x,q,y) = J [F{x + q,(p r x+ q y) - F(x)]dr ■ V q v = 

1-00 

= L [ F (x + q>fx+ q (<p S x+ q y))-F(x)]dr-V q v , 

J u 

where the group property of ^| is used in the second equality. Then, from 
the first identity in (|2.7[) , it follows that 



(2.8) —u{s,x,q,y) 



= - [F(x + q,y)~ F(x)] ■ V q v , 

s=0 



and from the second identity in (|2.7[) it follows that u(s, x, q, y) is in fact an 
explicit function of (pl +q y, that is, u(s,x,q,y) = U(x,q,(pl +(] y). However, 
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any U(x, q, tfz+qV) m ust obey the transport equation 

(2 9) ds 11 ^' q ' ^+1^ = VU (*' q ' ' = 

= G(x + q, <p x+q y) • VU(i, q, <p x+q y), 

where the second identity is due to (|2.6|) , and which holds for any s includ- 
ing s = 0. Combining (|2.8)> and (|2.9|) at s = 0, we obtain 

(2.10) G{x + q,y) ■ V y u{0,x,q,y) = -[F(x + q,y) -F(X)] ■ V q v . 

From the comparison with (|2.5b|) it follows that z?i = u(0, x, t/, t/), that is, 

/■OO 

(2.11) y 1= / [F{x + q f f x+q y)-F{x)]ds-VqV Q . 

J 

• Order e°. Here observe that vq does not depend on y, as pointed out above. 
This means that the average of Vo with respect to the invariant distribution 
measure }i x + q of (|1.3|) is the identity operation, and the same holds for its 
T-derivative. Then, averaging out (|2.5c|) with respect to }i x +q yields, 

lB = I n G{x + q,y)-V y v 2 d}i x+q {y) + 

+ J^ N}/ [F(x + q r y) - F(x)] ■ V q v r d}i x+q (y). 

For the first integral above we express the gradient term as a time deriva- 
tive of the function along the flow (in the same way as above for v\): 



(2.13) G(x + q,y) • V ' y v 2 (r,X,q,y) = —v 2 {T,x,q,(p% +q y) 



s=0 



The invariant distribution measure fi x +q preserves the averages of func- 
tions of f x+q y: 

(2.14) / v 2 (r,q,(p x+q y)d}i x+q (y) = constant, for all s. 
This, in turn, results in 

(2.15) — / v 2 (T,x,q,(p s x+CI y)d}i x+q (y) = 0, for all s, including s = 0, 

as Jr n j ' 

and, therefore, 

(2.16) ^ = J [F(x + q,y)- F(x)] ■ V q v x d}i x+q {y). 
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At this point, substituting the expression for x>\ from (|2.11|) , and rescaling time T 
back to t yields 



(2.17a) — 



Q(x + q,x)-V q + -S{x + q,X): ( V, <g> V,) 



vq, 



(2.17b) Q(a,b) = e J" (J-F(a,f a yfj [F{a,y) - F(b)]dji a (y)ds, 

[F(a,f a y)-F(b)] <g> [F(a,y) - F(b)]dp B (y)ds, 



(2.17c) 
S(a,b 



2eSym 



where ":" denotes the Frobenius product of two matrices, and "Sym" denotes 
the symmetric part of a matrix (the skew-symmetric part is canceled out by the 
Frobenius product with a symmetric matrix). 

As a result, the next-order correction q to the averaged dynamics in (|1.2|) , gen- 
erated by the backward Kolmogorov equation above, is given by the Ito diffusion 
process 

(2.18a) ^ = F(x), 

(2.18b) dq = Q(x + q, x)dt + cr(x + q, x)dW t , 

as long as q is small enough. Above, Wf is a K-dimensional Wiener process 
for some positive integer K, Q(x + q,x) is the 0(e) Ny-vector drift term, and 
<r(x + q,x) is the 0(^e) Ny x K stochastic diffusion matrix, given (non-uniquely) 
by 

(2.19) crcr T = S. 

Thus far, in (|2.18|) we obtained the equations for two sets of variables, x and q, 
while the evolution of the stochastic reduced dynamics for (|1.1[) is given by the 
sum (x + q). However, what we actually need is an approximate reduced model 
for (x + q) directly. For that, we add the two equations in (|2.18|) to obtain 

d(x + q) = F(x + q)dt + cr(x + q,x + q)dWt+ 

(2.20) +[F(x)-F(x + q) + Q(x + q, x)\ dt+ 

+ [cr(x + q,x) — cr(x + q,x + q)]dW t . 

Above, observe that the first term in the right-hand side is O(l), the second term 
is 0( y/e), and the rest of the terms are 0(e) or higher (taking into account that q is 
0(e)). Additionally, the terms in the second and third line of (|2.20[) are very diffi- 
cult to compute in practice for complex nonlinear dynamics. Therefore, we delete 
the terms in the second and third lines from (|2.20[) , thus completely decoupling 
(|2.20[) from (|2.18|) . Then, after replacing (x + q) — > x, we obtain the stochastic 
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reduced model in (|1.5|), where the diffusion matrix <r is computed according to 



(|2.17c|) and (|2.19|) with a = b = x (further we denote S(x) = S(x, x)). 



Remark: cr(x) does not depend on e. It is important to note that S(x), and, 
therefore, <x{x), do not in fact depend on e. Indeed, observe that 



(2.21) 



S(x) = 2eSym f f [F{x,<p x y) -F(x)] ® [F{x,y) - jF(*)]dft*(y)ds = 

= 2Sym f [F(x,<p s J E y) - F(x)] [F(x,y) - F(*)]d^(y)ds, 
Jo jR N y 

where (p s J £ y is the solution of the fast dynamics from (|1.1|) (with e _1 still in front of 
G(x, y)) with x = x, fixed as a constant parameter. Thus, there is no need to know 
what £ is to compute Six), and, in fact, £ does not have to be an explicit scaling 



parameter in the multiscale dynamics in (|1.1|) . This makes the new method prac- 
tical in realistic applications, where time-scale difference might not be explicitly 
available in the form of a parameter. 

3. Practical implementation of the reduced stochastic model for a 
general multiscale process with linear coupling 

As formulated above in Section |2l the new method is applicable for a broad 
range of dynamical systems with general forms of coupling. However, for the 
clarity of presentation, we dedicate a separate section to the implementation of 
the developed method for a multiscale process with linear coupling. The linear 
coupling is the most basic form of coupling in physical processes, however, be- 
cause of that it is also probably the most common form of coupling. Below we 
describe the complete step-by-step assembly of the reduced model, with both the 
deterministic and stochastic terms which parameterize coupling. 

Here we consider the special setting of (|1.1|) with linear coupling between x and 

y- 

dx dv 
(3-1) ^ = /(*) + L y y, -± = g (y) + LxX/ 

where / and g are nonlinear vector functions of x and y, respectively, and L x and 
Ly are constant matrices of suitable sizes. As before, x and y denote the slow 
and fast variables, respectively. However, one important distinction between (|3.1|) 
and (|1.1|) is that we no longer make use of the time scale separation parameter 
£; here the assumption is that it is somehow known that x-variables are slow, 
and t/-variables are fast, but no further information about time scale separation is 
available beyond that, and, in particular, no explicit time scale separation param- 
eter is known. As before, we introduce the corresponding fast limiting system 

dz 

(3-2) = g(z) + L x x 
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with x specified as a constant parameter. It is easy to see that for the multiscale 
system with linear coupling in (|3.1|) , the matrix S(x) from (|2.17c|) is computed as 
the time-lag correlation matrix 



(3.3) S(x) = 2LySym 



T 

RNy vPlv - *(*)) (y - *(*)) d^(y) ds 



where ^| is the solution of (|3.2[) , with x fixed at constant parameter x, forward 
in time s from an initial condition y, and z(x) is the mean state of (|3.2|) . Before 
constructing the reduced model, we need realistic assumptions on what we can 



compute in the multiscale system in (|3.1|) , which we take directly from HUH]]: 



First, we are going to presume that the multiscale dynamical system in (|3.1|) 
is not necessarily computable at will for arbitrarily long time intervals. The 
reason is that if it is possible, then the need for a reduced model becomes 
somewhat difficult to justify 

Even if the full multiscale model is not computable at will, we still need 
some statistical information about it to formulate the reduced model. Here, 
we presume that some typical state x* of the slow variables x is available, 
such that the dynamics evolve in the proximity of x*. For example, a 
rough estimate of the mean state of the slow variables of the full multiscale 
system can be taken as x* , or a nearby state. 



• We presume that the limiting fast dynamics in (|3.2[) is computable beyond 
the mixing time scale, so that time averages of (|3.2[) can be computed, at 
least for a single given value x = x* . 

• The assumptions above are the same as in BUS, in order to ensure com- 
patibility of the proposed stochastic reduced model method with what we 
have already developed in flUm for deterministic reduced models of mul- 
tiscale dynamics with nonlinear and multiplicative coupling. 

Under the above assumptions and the ergodicity hypothesis, we compute the 
mean state (z) and the time covariance matrix C(t) from the long-term time 
series of the solution z(f) of (|3.2|) with fixed parameter x = x*: 

(3.4a) (z) = lim - / z(t)dt, 

(3.4b) C(t)= Yirn\ ( T {z{t + T)-(z)){z{t)-(z)) T dt, 



(3.4c) C = / C(r)dT, 



CO 







Then, we assemble the reduced stochastic model for (|3.1|) in the explicit form as 



(3.5) dx = [f(x) + L y {z) + L y RL x (x - x*)] dt + crdW t , 
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Above, the deterministic part in (|3.5[) is computed as described in |H using the 
quasi-Gaussian linear response approximation ||27||, while the stochastic part is 
computed according to (|3.3|) (where the measure average is replaced with the time 
average) with x = x*. Observe that the matrix c will be computed only once, for 
the particular value x*. That, in effect, makes it a constant matrix approximation 
to the exact <x{x) from (|2.19|) , provided that the trajectory x(t) of the reduced 
model in (|3.5|) is in the vicinity of x*. This type of stochastic parameterization 
with a constant diffusion matrix is called the additive noise parameterization (as 
opposed to the multiplicative noise parameterization, where <r is a function of 
x). In the future work, we plan to extend the method described here onto the 
multiplicative noise parameterization. 

The choice of <x. Observe that c is not determined uniquely by (|3.6b|) , as multiple 
decompositions of S into crcr T are available (the Cholesky decomposition into the 
product of a lower-triangular matrix with its own transpose being one of the 
examples). Here we use the following algorithm to compute <r: first, solve the 
eigenvalue problem 



where X is the matrix of eigenvectors, and A is the diagonal matrix of eigenvalues. 
Since S is symmetric and nonnegative-definite II37L then all eigenvalues in A 
are guaranteed to be positive, and all eigenvectors in X are guaranteed to be 
orthonormal (so that X -1 = X T ). Then, we compute c as 



This method uniquely determines c as a square, symmetric, and positive-definite 
diffusion matrix for the Wiener process Wf. 

4. Computational study: the two-scale Lorenz 96 model with linear 



The test multiscale system used to study the new method of stochastic param- 
eterization here is the rescaled two-scale Lorenz 96 system with linear coupling, 
previously used in [4] to study the deterministic reduced model parameterization. 
The Lorenz 96 system is given by 



(3.7) 



SX = XA, 



(3.8) 



c = XA 1/2 X T . 



COUPLING 



(4.1a) ±{ = Xi_i(xi + i — x ( _2 



) - Xi) + 



F x -x 
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(4.1b) 

1 <.-.,,. s s , F y-y 



Vi,i 



y«,j+i(yy-i - y/,/+2) + i-{y{yi,j-i - y/j+2) - yy) + ^ 2 

Py Pi, 



£ 



with 1 < z < and 1 < / < /, such that Ny = N T /. The model has periodic 
boundary conditions: x !+ n v = x,-, yu+j = J/z+i,/7 and yi+N x ,j = yu- The parameter 
£ C 1 sets the time scale separation between the slow variables x ir and the fast 
variables y\u The parameters (x,p x ) and (y, j6y) are the (mean, standard devia- 
tion) pairs for the corresponding uncoupled and unrescaled Lorenz models 

(4.2a) xi = x f _i(x f+1 - x ! _ 2 ) - Xi + F x , 



(4-2b) yi t j = y i/j+1 (y i/hl - y iij+2 ) - y irj + F y/ 

with the same periodic boundary conditions. The rescaling above ensures that the 
Lorenz 96 model in (|4.1|) has zero mean state and unit standard deviation for both 
slow and fast variables in the absence of coupling (A T = Ay = 0), and remain near 
these values when X x and Ay are nonzero. The linear coupling above preserves 
the energy of the form 

\ N x f i N x } 

z i=l A J i=l 7=1 

Below we display the results of a numerical study of the stochastic reduced 
model for slow variables of multiscale dynamics, using the rescaled Lorenz sys- 



tem in (|4.1[) as the test model. We compare the statistical properties of the slow 
variables for the four following systems: 

(1) The complete rescaled Lorenz system from (|4.1|) ; 

(2) The stochastic reduced model from (|3.5[) ; 

(3) The deterministic reduced model obtained from (|3.5[) by removing the sto- 
chastic forcing (this model was previously developed in (H); 

(4) The poor man's version of (|3.5[) with no stochastic forcing and the first- 
order linear deterministic correction term R set to zero (further referred 
to as the "zero-order" system, same as in [4]). This zero-order system 
represents the simplest reduced model with constant parameterization of 
coupling terms. 

The fixed parameter x* for the computation of (z), R and tr was set to the long- 
term mean state (x) of the full multiscale model in (|4.1|) (in computationally ex- 
pensive models, a rough estimate could be used). Like in BHU, for all computa- 
tional statistical results presented below, the averaging time window T av equals 
10000 time units. 

Due to translational invariance of the studied models, the statistics are invariant 
with respect to the index shift for the variables X[. For diagnostics, we monitor 
the following long-term statistical quantities of Xf. 
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a. The distribution density functions, computed by bin-counting. A distribution 
density function gives the best information about the one-point statistics of X[, 
as it shows the statistical distribution of X{ in the phase space. 

b. The time auto-correlation functions (xj(f)xj(i + s)), where the time average is 
over t, normalized by the variance (xf) (so that it always starts with 1). 

c. The time cross-correlation functions (xi{t)xi + \{t + s)), also normalized by the 
variance (xj). 

d. The energy auto-correlation function 

w = (*?(Q*?(t + s)) 
1 ; (x2)2 + 2(x ! (i)^(i + s)) 2 ' 

This energy auto-correlation function measures the non-Gaussianity of the pro- 
cess (it is identically 1 for all s if the process is Gaussian, such as the Ornstein- 
Uhlenbeck process). For details, see [32 J. 

The following dynamical regimes are studied: 

• N x = 20, / = 4 (so that N v = 80). Thus, the number of the fast variables is 
four times greater than the number of the slow variables. 

• e = 0.01,0.1. The time scale separation of two orders of magnitude (e = 
0.01) is consistent with typical real-world geophysical processes (for ex- 
ample, the annual and diurnal cycles of the Earth's atmosphere). Also, in 
large-scale atmospheric processes the time scale separation between slow 
and fast variables can be weaker than that (for example, typical time scale 
of equatorial Kelvin waves is about 70 days, and that of Yanai waves is 
about 30 days, which is well in between the annual and diurnal cycles), so 
we additionally test the dynamical regimes with weak time scale separa- 
tion e = 0.1. 

• A x = A y = 0.3,0.35. These values of coupling are chosen so that they 
are neither too weak, nor too strong (although 0.3 is weaker, and 0.35 is 
stronger). However, this small variation in coupling changes the dynamical 
regime in the slow variables between moderately (e = 0.3) and weakly (e = 
0.35) chaotic and mixing, due to the suppression of chaos effect previously 
studied in 0. 

• F x = 6. The slow forcing F x adjusts the chaos and mixing properties of 
the slow variables, and in this work it is set to a weak-to-moderate chaotic 
regime F x = 6. The reason for that is that it was found previously in |4] that 
in strongly chaotic and mixing regimes at slow variables there is not much 
of a difference between the multiscale dynamics and reduced models. 

• Fy = 16. The fast forcing adjusts the chaos and mixing properties of the 
fast variables. Here the value of Fy is chosen so that the fast variables are 
strongly chaotic and mixing for Fy = 16. 
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Stochastic 


Deterministic 


Zero-order 


Density 


3.803 • 10" 3 


7.424 • 10" 3 


2.093 • 10" 2 


Corr. 


0.1218 


0.1152 


0.1935 


Cross-corr. 


0.1297 


0.1222 


0.2118 


Energy corr. 


1.312 • 10" 2 


1.436 • 10" 2 


3.473 • 10" 2 



Table 4.1. Relative errors between the slow variables of the full mul- 
tiscale system, and different reduced models, computed for the plots 
in Figure |4J] F x = 6, F y = 16, X x = A y = 0.3, e = 0.1. 



4.1. Moderate mixing at slow variables with weak time scale separation. Here 
we present the comparison of different statistics for the dynamical regime with 
moderate chaos and mixing at slow variables (achieved by setting A x = Ay = 0.3) 
and weak time scale separation (achieved by setting e = 0.1). In Figure 14.11 we 
show the distribution density functions, time auto-correlation functions, time 
cross-correlation functions, and energy auto-correlation functions for the multi- 
scale dynamics in (|4.1[) , and three different kinds of the reduced models: the 



new stochastic reduced model, the deterministic reduced model from |4j, and 
the zero-order reduced model with constant parameterization of coupling terms 
(which is given for reference as a simplest/poorest form of parameterization). 
Observe that the difference between the new stochastic and deterministic models 
here are rather small, however, it does look like the additional stochastic term 
improves the statistics. In particular, the stochastic terms makes the distribution 
density more spiky towards the multiscale density, and appears to produce bet- 
ter fits for correlation functions. The probable reason why there is not much of 
a difference between the deterministic and stochastic reduced models is that the 
slow dynamics are not sufficiently weakly chaotic for the stochastic term to make 
a significant difference. Table 14.11 confirms that there is some improvement error- 
wise in the distribution density and energy auto-correlation, but no improvement 
in auto- and cross-correlations (probably due to oscillations getting out of sync 
with increased lag time). 

4.2. Moderate mixing at slow variables with strong time scale separation. Here 
we present the comparison of different statistics for the dynamical regime with 
moderate chaos and mixing at slow variables (achieved by setting A T = Ay = 0.3) 
and strong time scale separation (achieved by setting e = 0.01). In Figure 14.21 
we show the distribution density functions, time auto-correlation functions, time 
cross-correlation functions, and energy auto-correlation functions for the multi- 
scale dynamics in (|4.1|), and three different kinds of the reduced models: the new 
stochastic reduced model, the deterministic reduced model from flU, and the zero- 
order reduced model with constant parameterization of coupling terms. Observe 



A SIMPLE STOCHASTIC PARAMETERIZATION FOR REDUCED MODELS 



15 



=20, N y =80, ^=0.3, Xy=0.3, F x =6, F y =16, e=0.1 



Q 



N x =20, N y =80, ^=0.3, ?i y =0.3, F x =6, F y =1 6, e=0.1 




=20, N y =80, ^=0.3, V - 3 ' F x= 6 ' F y= 16 ' e=0 - 1 
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Multiscale system 
Stochastic model 
Deterministic model 
Zero-order model 




Time 



Time 



Figure 4.1. Upper-left - distribution density function, upper-right 
- time auto-correlation correlation function, lower-left - time cross- 
correlation function, lower-right - energy auto-correlation function. 

F x = 6,F V = 16, X x = A y = 0.3, e = 0.1. ' 



that the difference between the new stochastic and deterministic models here are 
even smaller than for the regime with weak time-scale separation, probably due to 
the fact that the contribution of the stochastic term scales as y/e. However, again, 
it looks like the additional stochastic term improves the statistics somewhat. Table 
14.21 show that there is some improvement error-wise in all presented statistics, but 
not enough to be discernible visually in Figure 14.21 

4.3. Weak mixing at slow variables with weak time scale separation. Here we 
present the comparison of different statistics for the dynamical regime with weak 
chaos and mixing at slow variables (achieved by setting X x = A y = 0.35) and weak 
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Figure 4.2. Upper-left - distribution density function, upper-right 
- time auto-correlation correlation function, lower-left - time cross- 
correlation function, lower-right - energy auto-correlation function. 
F x = 6, F y = 16, X x = X y = 0.3, e = 0.01. 





Stochastic 


Deterministic 


Zero-order 


Density 
Corr. 
Cross-corr. 
Energy corr. 


8.105 • 10 -a 
9.309 • 10" 2 
9.57 • 10" 2 
9.042 • 10" 3 


1.048 • 10" z 
9.627 • 10" 2 
9.99 • 10" 2 
1.209 • 10" 2 


2.233 • 10~ 2 
0.1923 
0.2129 

2.776 • 10" 2 



Table 4.2. Relative errors between the slow variables of the full mul- 
tiscale system, and different reduced models, computed for the plots 
in Figure W^F x = 6,F y = 16, k x = X y = 0.3, e = 0.01. 
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N x =20, N v =80, 2u.=0.35, F x =6, F v =16, e=0.1 
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Figure 4.3. Upper-left - distribution density function, upper-right 
- time auto-correlation correlation function, lower-left - time cross- 
correlation function, lower-right - energy auto-correlation function. 
F x = 6,F V = 16, X x = A v = 0.35, e = 0.1. 





Stochastic 


Deterministic 


Zero-order 


Density 


2.166 • 10~ 2 


4.83 • 10~ 2 


7.516 • 10" 2 


Corr. 


0.2322 


0.2335 


0.3584 


Cross-corr. 


0.2277 


0.2346 


0.3557 


Energy corr. 


2.858 • 10" 2 


5.031 • 10- 2 


0.2163 



Table 4.3. Relative errors between the slow variables of the full mul- 
tiscale system, and different reduced models, computed for the plots 
in Figure EEH F x = 6, F y = 16, k x = A y = 0.35, e = 0.1. 
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Stochastic 


Deterministic 


Zero-order 


Density 


6.237 • 10" 2 


7.716 • 10" 2 


0.1088 


Corr. 


0.2629 


0.2752 


0.3769 


Cross-corr. 


0.2556 


0.2684 


0.3726 


Energy corr. 


0.1254 


0.1846 


0.3059 



Table 4.4. Relative errors between the slow variables of the full mul- 
tiscale system, and different reduced models, computed for the plots 
in Figure WM F x = 6, F y = 16, X x = \ y = 0.35, e = 0.01. 



time scale separation (achieved by setting e = 0.1). In Figure l4~3l we show the dis- 
tribution density functions, time auto-correlation functions, time cross-correlation 
functions, and energy auto-correlation functions for the multiscale dynamics in 
(|4.1|) , and three different kinds of the reduced models: the new stochastic reduced 
model, the deterministic reduced model from |H, and the zero-order reduced 
model with constant parameterization of coupling terms. Here we can see a sig- 
nificant improvement between the deterministic reduced model from [4] and the 
new stochastic reduced model. In particular, what apparently happens here is 
that the deterministic reduced model from [4j turns out to be less chaotic and 
mixing than the original multiscale dynamics at slow variables (observe that the 
distribution density is very spiky for the deterministic reduced model, while the 
time auto- and cross-correlation functions are less mixing, and the energy auto- 
correlation function is more sub-Gaussian than those for the multiscale dynam- 
ics). Then, the introduction of the stochastic term results in improvement of mix- 
ing and sub-Gaussianity, and also smoothens out the spikes on the distribution 
density, resulting in better approximation of the multiscale dynamics. Table 14.31 
confirms that there is improvement of statistics error-wise with the introduction 
of the stochastic term in the new reduced model. 

4.4. Weak mixing at slow variables with strong time scale separation. Here we 
present the comparison of different statistics for the dynamical regime with weak 
chaos and mixing at slow variables (achieved by setting A x = Ay = 0.35) and 
strong time scale separation (achieved by setting e = 0.01). In Figure 14.41 we 
show the distribution density functions, time auto-correlation functions, time 
cross-correlation functions, and energy auto-correlation functions for the multi- 
scale dynamics in (|4.1|) , and three different kinds of the reduced models: the new 
stochastic reduced model, the deterministic reduced model from ||4l, and the zero- 
order reduced model with constant parameterization of coupling terms. Here, 
again, we can see a significant improvement between the deterministic reduced 
model from [4] and the new stochastic reduced model, however, in a somewhat 
reverse way if compared to the previous set-up with weak time scale separation. 
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Figure 4.4. Upper-left - distribution density function, upper-right 
- time auto-correlation correlation function, lower-left - time cross- 
correlation function, lower-right - energy auto-correlation function. 



6, F v = 16, A 
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0.35, e = 0.01. 



Here observe that the deterministic reduced model is more chaotic and mixing 
than the slow variables of the multiscale dynamics, and the stochastic term makes 
the new model less chaotic and mixing to better match the multiscale dynamics. 
While at first seeming counter-intuitive, this effect can happen due to the sto- 
chastic noise pushing the solution off the unstable manifold of the deterministic 
system (where the chaotic and mixing motion occurs) into the dissipative absorb- 
ing region surrounding the system's attractor while not having enough strength 
to compensate for the lack of chaos and mixing with its own random forcing. In 
fact, the same (but somewhat weaker) effect can also be observed in Figures 14.11 
and 14.21 for weaker coupling and stronger chaos and mixing, so one can presume 
that it should be generally common in stochastic reduced models, especially if 



20 



RAFAIL V. ABRAMOV 



the time scale separation is strong and, because of that, the stochastic noise is 
weak enough to produce its own mixing. Table 14.41 confirms that there is some 
improvement of statistics error-wise with the introduction of the stochastic term 
in the new reduced model. 

5. Conclusions 

In the current work we improve the recently developed method HHU for de- 
terministic reduced models of multiscale dynamics with the higher-order addi- 
tive stochastic term which parameterizes the slow-fast interactions with random 
noise. The method is based on the homogenization techniques for multiscale sys- 
tems pTfl and offers a practical way of creating stochastic reduced models for a 
broad range of general multiscale processes. As demonstrated above in Section 
IH the stochastic term upgrade comes at no additional computational cost for sys- 
tems with linear coupling between slow and fast variables, as it uses the same 
time correlation matrix of the fast variables for a fixed state of the slow variables 
as the response term in the deterministic part of coupling parameterization. An- 
other practical advantage of the new method is that it does not require an explicit 
time-scale separation parameter between the slow and fast variables of multiscale 
dynamics. We tested the new method numerically using the two-scale Lorenz 
96 model with linear coupling between the slow and fast variables in a range 
of dynamical regimes with weak/ strong time-scale separation, chaos and mixing 
at the slow variables. The new stochastic reduced model consistently improved 
the results of the previously developed deterministic approach EHU in the two 
different dynamical regimes: 

(1) In the situation where the deterministic reduced model was less chaotic 
and weaker mixing than the slow variables of the full multiscale dynamics, 
the stochastic model was more chaotic and stronger mixing, due to stochas- 
tic forcing smoothening out spikes in distribution density and introducing 
random decorrelation in time auto- and cross-correlation functions. 

(2) In the situation where the deterministic reduced model was more chaotic 
and stronger mixing than the slow variables of the full multiscale dynam- 
ics, the stochastic model suppressed chaos and mixing in the reduced dy- 
namics. While this effect seems somewhat counter-intuitive, it can be ex- 
plained by the stochastic noise pushing the solution off the unstable man- 
ifold of the deterministic system (where the chaotic and mixing motion 
occurs) into the dissipative region around the system's attractor, while not 
having enough random force to increase chaos and mixing on its own. 

In the future work, we plan to extend the new method of stochastic parameter- 
ization of reduced slow dynamics onto multiplicative stochastic forcing param- 
eterization. Observe that the additive stochastic forcing parameterization in the 
current work emerges from the constant diffusion matrix approximation above in 
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Section |3] However, if the constant diffusion matrix approximation is improved by 
including its higher order Taylor expansion terms (as was done to the determin- 
istic coupling parameterization in IH[6]]), the inclusion of the higher-order terms 
will result in the multiplicative coupling with the diffusion matrix of the reduced 
model being a function of the slow variables. For the multiplicative diffusion 
matrix approximation, we plan to use an approach similar to that in dHS). 
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